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This paper investigates the energy harvested from the flutter of a plate in an axial flow by making 
use of piezoelectric materials. The equations for fully-coupled linear dynamics of the fluid-solid 
and electrical systems are derived. The continuous limit is then considered, when the characteristic 
length of the plate's deformations is large compared to the piezoelectric patches' length. The linear 
stability analysis of the coupled system is addressed from both a local and global point of view. 
Piezoelectric energy harvesting adds rigidity and damping on the motion of the flexible plate, and 
destabilization by dissipation is observed for negative energy waves propagating in the medium. 
This result is confirmed in the global analysis of fluttering modes of a finite-length plate. It is 
finally observed that waves or modes destabilized by piezoelectric coupling maximize the energy 
conversion efficiency. 



I. INTRODUCTION 

The environmental impact and limited resources of fossile fuel energies have motivated a significant research effort in 
the development of new and diverse techniques for the production of electrical energy. In parallel, a particular attention 
has also been given to systems able to produce limited amounts of energy at low cost to power remote or isolated 
devices, for which connection to the traditional electrical network is prohibitive in terms of cost or technical complexity 
[l[. These two elements have increased the attention on mechanisms able to produce self-sustained vibrations of a 
solid substrate on one hand and to convert the corresponding mechanical energy into electrical power on the other. 

The conversion into electricity of kinetic energy from geophysical flows such as tidal currents, winds and river 
flows is particularly attractive because of the large availability worldwide and the low environmental impact of this 
energy source 0. Research on fluid-solid interactions has identified several instability mechanisms that can lead 
to self-sustained vibrations of a solid placed in a steady uniform flow. Such fundamental instability mechanisms as 
coupled- mode flutter of a heaving and pitching airfoil vortex- induced vibrations Q or transverse galloping of 
flexibly mounted structures [j| are at the core of prototypes or concepts of flow energy harvesters. 

The harvesting of flow energy through flapping of thin clastic plates has also been investigated, mainly using two 
fundamentally different configurations. In the first one, an unsteady flow, created by the oscillatory wake of an 
upstream obstacle applies an unsteady forcing on the plate to make it flap 0, Q . The second configuration uses 
the coupled- mode flutter instability of the flexible plate in a steady flow to generate self-sustained flapping Q. In 
that case, it is well-known that the flat equilibrium state of the plate becomes unstable above a critical flow velocity, 
beyond which dynamic vibrations of large amplitude can develop on the structure. The linear stability of this system 
has been extensively considered from both a local and a global point of view. In the former, the instability of waves 
in the infinite medium was considered ffjl. [HI , w hile the latter considered finite systems, including effects such as 
vortex shedding downstream of the plate [ll[ [l2j , three-dimensional effects [l3| , lateral confinement [l4| , spanwise 
confinement [la ] or coupling between multiple structures [l6j . The non-linear self-sustained flapping devel oping above 
the instability threshold has also been the focus of multiple experimental [lT], EH and numerical studies (l9rl2l1 ] . 

In this work, we are interested in the ability to produce electrical power from the self-sustained oscillations of a 
flexible plate resulting from this fluttering instability. To assess the potential for electrical energy production, it 
is important to properly include in the dynamical equations of the fluid-solid sytem the loss of energy due to the 
conversion into electricity, as we seek here regimes where the extracted energy is a significant fraction of the total 
energy of the system. Two main approaches can be considered to produce electricity from the mechanical energy 
of a vibrating solid. Classical generators convert a displacement of a solid substrate into electrical energy through 
electromagnetic induction, and are commonly used in classical turbines as well as in recent prototypes of flow energy 
harvesters On the other hand, piezo-electric materials convert mechanical strain into electric potential, and have 
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recently received increasing attention for applications involving low power production, typically of the order of the 
mW [l[ HH . Studies on piezoelectric materials for energy harvesting considered flow- induced vibrations [f| H^, HU but 
also vibrations from various other forcings such as human movements pH . [26| • 

The objective of the present work is to study from a theoretical point of view the stability properties and dynamics of 
a classical fluid-solid system (a flexible plate subject to coupled-mode flutter) coupled to an output electrical network 
with piezoelectric materials. In comparison to previous studies on flow energy harvesting Q or piezoelectric energy 
conversion, the present approach is original by its full coupling of the fluid-solid and electrical systems. In this paper, 
the linear stability of the fluid-solid-clectric system is investigated from both a local and global point of view to identify 
the impact of the coupling on the stability of the system and assess the efficiency of the mechanical-to-electrical energy 
conversion. 

Coupling the vibrations with electrical circuits that dissipate energy intuitively results in damping from the point 
of view of the structure. In the context of fluid-structure interactions, the effect of viscous or viscoelastic damping has 
been addressed both on the infinite length medium (local approach) and the finite length systems (global approach). 
[lol ] investigated the effect of damping on the stability of flexural waves propagating in compliant panels interacting 
with an homogeneous flow and found that dissipation can destabilize some particular waves, which are referred to 
as Negative Energy Waves after [27| . In the finite length case, the most studied system with respect to the effect of 
damping is the fluid- conveying pipe, which shares many similarities with the plate in axial flow (28l. |29|. The effect 
of p iezoelectric coupling on the instability threshold of a cantilevered fluid- conveying pipe has been addressed by 
[301 ] . Destabilization or stabilization by dissipation has been observed, depending on the fluid-solid mass-ratio. The 
comparison of local and global instability criteria with damping has been addressed by 31| . It was shown that global 
instability of the long system is always predicted by the local instability criterion of the dissipative medium. 

The particular effect of damping induced by piezoelectric coupling has been widely investigated in the research field 
of structural damping and vibration control. Passive damping by the use of shunted piezoelectric patches, developed 
by 1321] . has inspired studies involving piezoelectric materials and passive electrical components arranged in a network 
|33l |34| . We will consider here the continuous limit addressed in these two studies, that is valid when the length of 
the piezoelectric material is small compared to the typical wavelength of the solid's deformations. 

In Section UH the linearized equations of motion for the coupled piezo- mechanical problem of a flexible plate covered 
with piezo-electric elements in a potential flow will be presented. In Sections IIIII and IIV[ the linear stability analysis 
will be investigated from a local and global point of view, respectively. In both cases, the effect of piezoelectric 
coupling on stability will be investigated, as well as the efficiency of energy conversion between the solid mechanical 
energy and the energy dissipated in a simple resistive electrical network. Finally, in Section [VI the local and global 
stability results will be discussed as well as possible extensions of the present work. 



II. PROBLEM FORMULATION 



We consider here the motion of a flexible plate of thickness ho, Young's modulus Eq, density po and Poisson's coefficient 
vq- The span of the flexible plate is assumed to be much greater than its typical streamwise lengthscale. We therefore 
focus on purely two-dimensional deformations, and the plate's vertical displacement is noted w(x,t). Piezoelectric 
patches of length I, thickness h p , density p p , electric permeability e, Young's modulus E p and Poisson's coefficient v p 
are attached symmetrically on the plate. The plate is surrounded by a fluid of density pf with upstream horizontal 
velocity Uoo . The problem is sketched in Figure [T] In this section, we present the linearized coupled equations of the 
fluid-solid and piezo-electric systems, first in the case of discrete piezoelectric patches, and then in the continuous 
limit when I is much smaller than the streamwise lengthscale of the solid deformation. In this limit, the equation 
of energy conservation will enable us to exhibit the different energy transfers in the system. In the following, for a 
function a(x,t), a and a' correspond to the temporal and streamwise derivatives of a respectively. 



A. Equilibrium equations of a beam with discrete pairs of piezoelectric patches 

The left and right ends of the i piezoelectric pair are positioned in xj and xf, so that I = xf — x~ and 
Xi = (x~[ + xf)/2 denotes the position of the center of the patch. We consider complete coverage of the plate by 
piezoelectric patches, therefore xf = x~[, x . The following equations can easily be adapted to the case of partial 
coverage [34[ . Quantities related to the piezoelectric patch located on the upper (respectively, lower) face are denoted 
by the exponent W (respectively, ^). The patches are attached to the plate so that their respective polarities are 
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FIG. 1: Schematic view of a plate in an homogeneous axial flow, equipped with small length piezoelectric patches on both 
sides. 

reversed, and the charge displacement across a piezoelectric patch is obtained as [H, H(| : 



,(fe) 



X 



Cv} k \ k = l,2, 



(1) 



where \ = e3i(/io + /i p )/2 is a mechanical/piezoelectrical conversion factor with e^i the coupling factor, and C = el/h p 

is the capacity per unit length in the spanwise direction of the piezoelectric element. is the potential difference 
between the two electrodes of the corresponding piezoelectric (see Figure [T]) . Negative electrodes of the patches are 
assumed to be shunted through the plate by a purely conductive material, therefore, 



and 
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V,.. 



(2) 



(3) 



Each piezoelectric pair is therefore equivalent to a single piezoelectric patch of equivalent capacity C and voltage V% 
of respective values, 



C=j 7 V q =2V l7 



(4) 



and this representation is retained for simplicity in the remaining of the paper. Assuming an Eulcr-Bcrnoulli model 
for the plate, the bending moment at a given position x results from the internal rigidity of the material and the 
piezoelectric coupling: 



M = Bw" -J2xV t F t . 



(5) 



where B is the flexural rigidity of the three-layer plate [3 

12(1^ 2 ) + 1 
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and Fi is the characteristic function of the i th patch 

Fi = H(x-xr)-H(x-xf). 



(7) 
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FIG. 2: All piezoelectric pairs are shunted with resistances, modeling the electrical energy absorption. 



with H the Heavyside step function. In Eq. (JSJ, the summation over all piezoelectric patches and the use of the 
characteristic functions Fi provides a compact expression for the local bending moment induced by the piezoelectric 
coupling. The linearized conservation of momentum for the plate then leads to 



-M"-[P], 



(8) 



where \i is the surface density of the plate with piezoelectric elements and [P] = P + — P~ is the pressure forcing 
applied by the fluid on the plate. We consider here a purely inviscid model for the flow, and the exact expression of [P] 
will be made explicit in the following sections for the infinite- and finite-length systems respectively. In non-inviscid 
flow, the viscous boundary layers would tend to stabilize the plate by inducing a non-uniform tension, maximum 
upstream [see for example 1 2 lj. For large enough Reynolds numbers, this correction is however expected to be small. 

An additional relation between Vj and Qi is needed to close the system of equations (JTJ, © and (El) , and is obtained 
from the output electrical network connected to the free electrodes of the piezoelectric patches. [38[ showed that 
energy removal from the piezoelectric system induces structural dissipation, whatever electric device is effectively used 
(resistor, storage in a battery or other energy harvesting circuitry). Consequently, the simplest electrical component, 
where each piezoelectric pair is shunted by a conductance G, is considered to investigate the potential harvesting 
power of the system and its effect on the flutter dynamics (Figure [5]). Applying Ohm's law, 



Vi 



Qt 
G 



= 0. 



This simple electric network is a particular case of more general electric networks considered in [3J 
the optimization of dynamical properties of phmcd-phmcd beams. 



(9) 



in the context of 



B. Continuous limit 



If the typical lengthscale of the plate's streamwise deformations is large compared to the length I of the piezoelectric 
patches, one may consider the continuous limit of the discrete equations JTJ), ©, © an d ©■ We introduce the surface 
density of the charge, piezoelectric capacity and conductance of the resistive circuit as 



q{xi) = Qi/l, c=C/l, g = G/l, 
and the continuous voltage v(x = Xi) = Vi- As I — > 0, 



["'' 



w"(Xi)l, 

v(x). 



(10) 

(11) 
(12) 



The continuous equations are then obtained as 



/iw + Bw"" — x v " = — [P] ; 
q — cv — \w" = 0, 
F(q,v) = 0. 



(13) 
(14) 
(15) 



Here, F(q, v) is a general expression that links charge and voltage, so that at this stage, the equations are valid for 
any circuit. In the particular case considered here, F(q, v) = v + q/g and v can be eliminated from the previous 
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dynamical equations, thereby leading to a system for w and q only: 
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-q + q — X w = 0. 
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(16) 
(17) 



These continuous equations are similar to those obtained in the Laplace space by [34| where homogenization techniques 
were used. 



C. Efficiency of the energy conversion 



To assess the amount of electrical power that can effectively be extracted, nonlinear effects are important to provide 
the saturation amplitude of the self-sustained oscillations. This question is not addressed here as we focus only on 
the linear analysis, and will be the subject of a subsequent contribution. However, important physical insight can 
be gained from the analysis of energy transfer by linearly unstable modes. In particular, the non-linear mode shape 
that determines the piezoelectric deformation rate and the e nerg y transfers to the output circuit, has been observed 
to be similar to that of the most linearly unstable mode [2fj,|39(. From Eqs. (fT5)l and (|14[) . the equation for energy 
conservation can be obtained as 
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with 

l B w" 2 + ±cv 2 , (19) 

Bit/") + w' (Bw" - xv) , (20) 

(21) 
(22) 

The total energy density of the system $ is the sum of the solid kinetic and elastic energy as well as the electrical 
energy stored in the capacity of the piezoelectric material. & is the flux of mechanical energy in the plate: the first 
and second terms are respectively the rate of work of the internal bending forces and moments, due both to the elastic 
response of the solid and the piezoelectric coupling. 3P p is the rate of work of the local pressure force and SPei is the 
power readily available and dissipated in the output system. Note that this equation is valid regardless of the output 
electrical network chosen. In the case of the purely resistive system considered here, = q 2 / 1 g. 

We consider a measure of the harvested (i.e. dissipated) energy over one flapping period T relative to the mean 
energy density of the solid-piezoelectric system during that period. In that regard, this ratio is a normalization of 
the harvested energy, and it will be referred to as conversion efficiency in the following. This ratio is defined for an 
unstable mode (that can lead to self-sustained oscillations) as 

T 



o 



, ,T ■ (23) 

? / 



In this last expression (.) stands for the spatial mean value for the considered mode, taken over either a wavelength 
in the local analysis or the entire plate in the global analysis. Note that since r is just a normalized energy output 
and not a thermodynamic efficiency, r > 1 is allowed. In the following, we will study the influence of the system's 
parameters on this ratio to find optimal conditions for the energy conversion. 
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III. DESTABILIZATION BY DAMPING AND ENERGY CONVERSION IN THE INFINITE MEDIUM 

A. Non-dimensional equations 

Using p/pf, fi/pfUoo, pU^a and Uooy/pc respectively as characteristic length, time, pressure and charge surface density, 
and noting all non-dimensional variables with a tilde, Eqs. (|16p and (|17p take the following non-dimensional form: 
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(24) 
(25) 

(26) 
(27) 
(28) 



V* is the non-dimensional flow velocity, a is the piezoelectric coupling coefficient, and 7 is the ratio between the 
fluid-solid and electrical characteristic timescales. 

In the present local analysis, mechanical and electrical displacements are sought in the form of harmonic waves of 
wavenumber k and frequency a), 

(29) 



a; 


= Re 








A 1o J 



B. Computation of the pressure forcing 



The pressure forces are computed assuming a potential flow on both sides of the flexible solid. The potential $ can 
be decomposed in ^(a;, y, t) — x + 4>{x, y, t) where </> is the perturbation potential to the uniform base flow. The flow 
is incompressible therefore <j) must satisfy 



with boundary conditions, 



V 2 = O, 



— {x,y = Q^,t) = —{x,y = Q ) = w{x, t) + w (x, t), 
-> 0, for \y\ -> 00. 



(30) 

(31) 
(32) 



Note that <j) is discontinuous on the plate but the normal velocity must be continuous. The pressure is obtained from 
the linearized unsteady Bernoulli equation as 



d(j) d(j) 
^ dt dx 



(33) 



The velocity potential is obtained solving Eq. (|3TJ)) with boundary conditions (|3"Tj) - (|3"2"|) . The pressure jump is then 
obtained by applying Eq. (f33|) on both sides of the plate. From Eq. (|29|) , it can be written 



2w 



(Q - k) 2 i(kx-ut) 



[p}=p(x,y = + ,t) -p(x,y = 0-,t)=Be 
Using Eqs. (|34|) and (f29| . the system (|24|) -(|25 |) becomes a linear system for wo and qo 

= C(k,Cj) 
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with 

D »=- i2 ( i+ it)- 2 « +4 4 + ^ (3e) 

D 2 = 1 - Z7W. (38) 

Equation (1351) admits non trivial solutions if and only if the determinant of C vanishes, which leads to the following 
dispersion relation: 

D (k,uj) = D 1 (k,u>f(— -lV (39) 

For (k, ui) solution of Eq. (|39f) . Eq. (|35[) determines the relative amplitude of the mechanical displacement and electrical 
charge for the corresponding wave. 



C. Local stability analysis 



We will restrict the wave analysis to the temporal approach, which classically consists in looking for frequencies 
satisfying the dispersion relation associated to a real wavenumber k. Equation Q39[) can be put in the form of a 
third-order polynomial in uj. Hence, for any real wavenumber, there are three different waves. If there exists a 
real wavenumber for which one of these frequencies has a positive imaginary part, the corresponding wave grows 
exponentially in time, indicating a temporal instability. In the energy harvesting context, this temporal instability 
is necessary to create self-sustained oscillations of the plate that are able to generate a net power to the electric 
networks via the piezoelectric coupling. 

Let us first consider the case of no-piezoelectric coupling (a = 0): the matrix C is diagonal and the problem then 
consists in two distinct dispersion relations. The first describes the propagation of flexural waves in the medium, 

D (k,oj) = 0. (40) 

This dispersion relation is very similar to that of a compliant panel interacting with a potential flow, which has been 
extensively studied 0, |T(| ■ The difference between the present dispersion relation and that of Q comes from the 
different choice for the non-dimensional time and the presence of the flow on both sides of the plate. Hence, the same 
phenomena will be observed, but at different values of the non-dimensional velocity, frequencies and wavenumbers. 
The main feature of this medium is that it is unstable for any non-zero value of the flow velocity. Analyses of 
the different branches in the complex k- and ^-planes have also been performed in the above-mentioned papers to 
investigate the convective-absolute instability transition, which will not be addressed here. 

The second uncoupled dispersion relation describes the behavior of the charge in the electric network, 

D 2 (k,u) = 0. (41) 

Equation (|4ip corresponds to the charge dynamics in an RC circuit and does not include the wavenumber k, as no 
propagation of charge can exist in a medium composed of electrical circuits disconnected from each other. 

On Figure[3l the frequencies Cj n {k) (n = 1..3) are plotted. Here, G>\ and u>2 are solutions of Eq. (j40|) and correspond 
to flexural waves, while cl>3 is the solution of Eq. (|4ip and corresponds to an electrical wave. For k £ [0, fc&], the two 
frequencies &i and (I>2 are complex conjugate, one of them having a positive imaginary part that indicates a temporal 
instability. For k > fej,, uj\ and W2 are real and the waves are neutrally stable. Wave 1 has positive phase velocity for 
all k but two ranges of wavenumbers can be distinguished depending on the sign of the phase velocity of wave 2. For 
kb < k < k c , LJ2 > and wave 2 has a positive phase velocity, while for k > k c , it has a negative phase velocity. As a 
consequence of the phase velocity sign change, the frequency Co vanishes for k = k c . From Eq. (|40p . one can show that 

h~V*, (42) 

when < 1. Looking for zeroes of u>2, we obtain 



(43) 



8 



3_ 





(c) 



- 

-0.01 - 

-0.02 - 

£ -0.03 - 
£ 

-0.04 - 

-0.05 - 

-0.06 - 

-0.07 ~ 




0.15 
k 



FIG. 3: Real part (a) and imaginary part (b,c) of frequencies associated to a real wavenumber k, for V* — 0.05, 7 = 15 and 
a — 0. This corresponds to a situation where no coupling is present in the system. 



Finally, the frequency C03 is constant and purely imaginary with negative imaginary part: no energy propagates and 
the negative growth rate Imijli^) is the characteristic time of discharge of a capacity c in a resistance l/g. 

Piezoelectric coupling is now considered (a / 0). From a mechanical point of view, its effect is to pump energy from 
the system to produce heat in the resistances, thereby dissipating mechanical energy. It is thus expected that results 
regarding the effect of damping on stability of neutral waves will also apply here. In particular, as demonstrated by 
[271 ] , the sign of wave energy is expected to allow the prediction of the stabilizing or destabilizing effect of coupling 
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FIG. 4: Wave energy of the neutral waves propagating in the system without piezoelectric coupling for V* = 0.05. 



on neutral waves. The wave's energy is defined as the work done in slowly building up the wave starting from rest at 
time t = — oo and has for expression [io| . 



= ui dDp 
4 du 



w 1 



(44) 



Let us consider a neutral wave propagating in the uncoupled system. The corresponding frequency is u> (k), so that 
Z?o(^Oj k) = and Co^i = 0. The perturbed value of Q due to the addition of piezoelectric coupling is then written as, 



LJ = LUq + Sill, 

where SCj <C Co. The frequency uj satisfies the dispersion relation, thus, 

dD 



SCj 



dCj 



(45) 



(46) 



(k,ui ) 



In this last expression, only the leading order terms have been kept. The stabilizing or destabilizing effect of piezo- 
electric coupling depends on the imaginary part of 5Cj. After a straightforward calculation, this quantity reads, 



Suji 



V* 2 (l + wg7 2 ) 



(47) 



It hence appears that the variation of the growth rate after the addition of piezoelectric coupling has the opposite sign 
of the neutral waves energy in absence of coupling. It is positive for a negative energy wave (NEW), and negative for 
a positive energy wave (PEW). Following the classification of [4l|, NEW are also referred to as class A waves, while 
PEW are referred to as class B waves. On FigureHl energy of waves 1 and 2 are plotted as function of the wavenumbcr. 
Energy of wave 2 is negative for k € [fc&, k c ] and the wave energy analysis shows that this range of wavenumbers will 
become unstable when piezoelectric coupling is added. To address the validity of the previous prediction, Figure [5] 
represents the three frequencies associated with a real wavenumbcr in the case V* = 0.05, a = 0.5, so that it is the 
same case as in Figures |3] and 31 but with added piezoelectric damping. Here, piezoelectricity couples mechanical 
displacement and electrical displacement. Consequently, no purely mechanical or electrical waves propagate in the 
system. The predictions of the above wave energy analysis are confirmed by the behavior of Im(a)i) and Im(d)2) on 
Figure [SJa: Wave 1 is stabilized by the addition of piezoelectric coupling while wave 2 is destabilized by the addition 
of piezoelectric coupling in the range of wavenumbers [kb,k c ] and stabilized for k > k c . For the third wave, which 
is a purely electrical wave without coupling, the frequency has now a non-zero real part, indicating that it is now a 
propagating wave. Consequently, this wave also appears to be affected by coupling. 
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FIG. 5: Real part (a) and imaginary part (b,c) of frequencies associated to a real wavenumber k, for V* — 0.05, 7 = 15 and 
a — 0.5. This corresponds to a situation where coupling is present in the system. 

D. Energy conversion efficiency of unstable waves 

Now that the effect of piezoelectric coupling on the propagation of waves and their stability has been addressed, 
the energy conversion efficiency of these waves is investigated. As defined in the previous section, energy conversion 
is modeled by resistances that shunts the piezoelectric patches, and is significant only if a wave is unstable, so 
that it grows exponentially in time and eventually saturates through a nonlinear mechanism that is not addressed 
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in the present analysis. The starting point is the ratio defined in Eq. f[23|) . which takes the following form when 
non-dimensional variables are used: 



1 

7 Jo 



2n/u> r 

(vq) dt 



~ />27r/tJ r 1 i 

tJ ^ 2 + ^" 2+i2 > dt 



(48) 



The spatial averaging appearing in this last expression corresponds to the mean value over one wavelength. When 
it is applied to the product of two quantities transported by a wave of wavenumber k and frequency u, namely 
a = Re[a e 1 ^-"*)] an d b = Refoe^ - "*)], it reads 

(ab) = — abdx = 2e 2 ^ t Rc(a b ). (49) 

2n Jo 

The numerator and denominator of Eq. (|48p have the same time-dependence therefore r does not depend on time. 
With a displacement in the form of Eq. (|29fl , it finally takes the following form, 



r(a,V*,-y,k,uj) = ^ 




V 



2 1 - 1 

?7W — 1 



vA -2 



+ 7 



(50) 



where Q and k are linked through the relation dispersion Eq. For a given set of parameters a, V* and 7, let R 

be the maximum value of r among all unstable waves, 

R(a, V* , 7) = max r(a, V*, 7, k, £>). (51) 

k.uii>0 

The wavenumber and frequency of the corresponding wave are noted K and Q, respectively. The maximum efficiency 
R is plotted in Fig. [5] as function of 7, for a = 0.5 and V* = 0.05. It presents a local maximum at 7 ~ 22 and tends 
to infinity for high values of 7. The values of K and O r are also plotted. The first observation that can be done is 
that for all 7, K is comprised between k^ and k c . Hence, the maximum of r always occurs for a wave that is stable 
without piezoelectric damping and is destabilized by addition of piezoelectric damping. 

Let us focus at first on the behavior of R for 7 — > 00. One is tempted to conclude that one has to choose a large 
value of 7 to optimize the energy conversion. But it appears on Figs. |6|b)-(c) that K — > k c and fl r — > u>2(k c ) = 0. 
This means that although the energy dissipated in the electrical circuits grows with 7, the period of the oscillations 
diverges and the growth rate tends to zero. This situation is thus far from optimal. One can also note that 7^00 
corresponds to a resistance that tends to infinity, or equivalently, an open circuit. For this reason, we conclude that 
the optimal efficiency corresponds to the local maximum of Figure |6^i around 7 = 22 instead of 7 — > 00. Moreover, 
Figure ISJc) shows that in the vicinity of the maximum efficiency, the characteristic time of the electrical circuits and 
the characteristic time of the wave are equal, indicating that optimal efficiency also results from a synchronization 
between the fluid-solid and electrical systems. 

For the particular values of the parameters a and V* used in Fig. |6l negative energy waves are hence waves that 
optimize conversion efficiency. Let us now explore the whole (a, V*) space. To do so, we introduce r(V*, a), the value 
of 7 that maximizes R(a,V* ,7) and i? 7 , K n and f2 7 the respective values of the efficiency, the wavenumber and the 
frequency at this maximum. On Fig. [71 K 1 is plotted as function V* for several values of a between 0.01 and 0.7. 
The range of wavenumbers destabilized by damping, [kb 7 k c ] appears grayed out on this plot, showing that for any 
values of a and V*, optimal efficiency occurs for a wave destabilized by damping. 

Finally, the optimal efficiency i? 7 is plotted as function of a for different values of V* on Fig. [5] All curves gather 
on a single line, indicating that the optimal efficiency scales as a and does not depend on the flow velocity. This 
enlightens the importance of maximizing the coupling coefficient for such application. 

Before addressing the finite length system, let us summarize the main point of the local analysis conducted in this 
section: The particular range k c ] of wavenumbers has been emphasized, in which piezoelectric damping has a 
destabilizing effect and analysis of the efficiency has shown that energy conversion efficiency is maximum for waves 
in this range. 



IV. LINEAR DYNAMICS OF A FINITE-LENGTH PIEZOELECTRIC FLEXIBLE PLATE 



In the previous section, we focused on the local stability analysis of an infinitely long flexible plate coupled through 
piezo-electric patches to a purely dissipative electrical system. In this section, we are interested in the global stability 
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FIG. 6: (a) Maximal efficiency, R as function of 7 for a = 0.5 and V* = 0.05; (b) Value of the wavenumber of the corresponding 
unstable wave; (c) value of the corresponding frequency, compared with I/7 (dashed line). 

analysis and conversion efficiency of a plate with finite length, coupled to the same piezo-electric dissipative system. 
The wave analysis approach of the local study is here replaced by a study of global modes that takes into account the 
boundary conditions on the finite plate and effects such as vortex shedding downstream from the solid. 
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FIG. 7: Value of the wavenumber that maximize the efficiency, K-y. Grayed region corresponds to the region of negative energy 
waves, k £ [kb,k c ], showing that these waves optimize the efficiency. 
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FIG. 8: Optimal energy conversion efficiency, i? 7 , as function of a for different values of the non-dimensional velocity V* . 

A. Non-dimensional equations and fluid model 

The natural length scale of the problem is now L, the streamwise length of the plate. The non-dimensionalization 
of the system's equations must be modified, and L and L/Uoo are now used as characteristic length and time scales 
respectively. Pressure, charge densities and voltage remain non-dimensionalized by pU^, U^^/JIc and U oc y /r jI/c. 
The non-dimensional form of a variable a is noted a to avoid confusion with the infinite system's notations, so that 
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Eqs. (lT6l) - (fT7|) now become: 



JJ* 



: {l + a z )w"" + w 



U* 



-M*p, 



where the coupling parameter a was defined in Eq. (|27[) and 



M* 



PfL 



7 

M* 



U* = UL\ 



V*M* 



The voltage induced on the piezoelectric system is in non-dimensional form: 

a „„ 
v = q-—w . 



(52) 
(53) 

(54) 
(55) 



The parameters M* and U* are the classical mass ratio and non-dimensional velocity used in previous flag stability 
analysis [HG1. P 

is the non-dimensional characteristic time-scale of the output circuit, relative to the advective 
time-scale L/Uoo chosen as a reference in this section. 

M* can also be interpreted as a length ratio L/r] where r\ = fi/pf is the characteristic length scale introduced in 
Section inTl The parameters U* and (3 are the finite-system equivalent to the V* and 7 parameters of the local analysis 
[see Eqs. (f2"6")l and (|28p]. In the limit of a long flexible plate, M* ^> 1, the local dynamics studied in the previous section 
is expected to become dominant over finite-length effects such as the influence of the wake or boundary conditions. 

Clamped-free boundary conditions arc imposed on the finite-length flexible plate: 



for x = 
for x = 1 



w 
w' 





= 



(1 + a 2 )w" - aU*q = 
(l + a 2 )w"' ~aU*q' = ' 



(56) 
(57) 



Note that, Eq. (f5"7} expresses the vanishing at the free end of the total internal torque and normal stress. 
The solution of the linearized equations Eqs. (f5"2|) - ([53|) is sought in the form of global modes 



Re 



/ W{x) 



(58) 



where o> can be complex. Similarly, V(x) can be defined from the potential v and is easily obtained from Eq. (|55[) . 



B. Pressure forcing on the finite- length plate 



The pressure forcing on the flexible plate is computed from a potential flow approximation. In the case of a finite- 
length solid, we must account for the shedding of a vortex wake from the trailing edge of the plate. We follow here 
the double- wake method [l3l. Il4l fltjj : The flow around the plate is potential everywhere except on the horizontal axis. 
On the plate, a bound vorticity distribution is present to satisfy the no- normal flow boundary condition on the plate's 
upper and lower surfaces. The free wake vorticity is not solved for explicitely as in the Vortex Sheet approach [TJ, [42[ 
but it is instead imposed that the pressure discontinuity across the plate must vanish at both ends of the flexible plate 
([p](0, = [p](l,t) = 0). At the trailing edge, this is consistent with the shedding of a free horizontal vortex sheet 
that can not sustain any pressure force. At the leading edge, it introduces an "upstream wake" whose physical origin 
is debatable. However, this method has been shown to predict the instability threshold and growth-rates correctly, 
particularly for intermediate-to- large mass ratio M* where the instability threshold and mode structures predicted 
by both methods are very similar [lfj . 

Using Fourier decomposition in the axial direction, it can be shown that the longitudinal gradient of the pressure 
jump d\p\/dx satisfies the following singular integral equation [l3j : 

SjfT?«'"i^-(! + s)'*<** H<M>-MM-* (59) 

This singular equation can be formally solved for d[p]/dx, and after integration of the pressure gradient from the 
leading edge, the pressure forcing on the plate is obtained as a function of the plate's displacement w(x, t). Using the 
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normal mode decomposition in Eq. (|58p. for a given mode shape W(x), the pressure jump P(x) can be decomposed 
as 

P{x) = -lu 2 P ai (x) - 2iuP G {x) + P K {x) (60) 

with 

r ig(P K m\ At [W"(x)\ ( P K [W]\ 

o~ WA PG[W] \ t=\ W ' {x) > PG [^ (» = 0,1) = 0. (61) 

P M , P G and P K represent added inertia, gyroscopic and added stiffness effects on the particular mode considered, 
and can be expressed as linear operators on the mode shape W by inverting Eq. (|6"Tj). The system (|52 p ~(|53 ]l can then 
be rewritten as 



-C? (W + P M [W]) - 2iuP G [W] + { ^" ; {W"" +P K [W}) - —Q" = (62) 

(-itS/3 + l)Q - jj-W" = 0, (63) 

and boundary conditions are readily obtained from Eqs. (|56[) - (|57|) . 

For given values of the four non-dimensional parameters a, /3, U* and M* , discretizing W(x) and Q(x) on the hrst 
N Chebyshev Gauss-Lobatto points, the above system together with boundary conditions in Eq. (|55 )) - ([57|) . can be 
written using a collocation method as an eigenvalue problem for the vector [W, — itliW, Q] T . 



C. Stability analysis and impact of piezo-electric coupling 

The flexible plate in axial flow classically becomes unstable to fluttering above a critical velocity ratio U* lit that 
depends on the fluid-solid mass ratio M* [e.g. [13]]. We consider here the evolution of this stability threshold with 
the piczo-clectric coupling a and the circuit's response time-scale /3. 

For a fixed (3, increasing a enhances the coupling between the piezo-electric and mechanical systems: an additional 
rigidity is introduced and energy is dissipated in the resistive circuit. The overall effect of the piezo-electric coupling 
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a is therefore expected to be a stabilization of the fluttering modes. Figure [9] shows that this is generally the case 
except at large values of M* where destabilization of some modes through the piezo-electric coupling is possible. For 
large enough a however, all modes are stabilized, in agreement with physical intuition on the stabilizing effect of a 
dissipative system. 

Nonetheless, it is important to note that the effect of a is highly dependent on the tuning of the system's frequency 
to that of the output circuit and Figure [9] corresponds to a configuration where both the fluid-solid system and the 
electrical circuit have similar fundamental frequencies. 

The ratio /3 is also a measure of the resistance of the electric circuit, and the limit /3 <C 1 corresponds to shunted 
piezo-electric elements. In that case, charge transfers are instantaneous and the electrostatic balance of the piezo- 
electric material is achieved at all times. The electric potential difference v in the piezo-electric remains negligible, and 
so does the piezo-electric torque on the plate. Therefore, in the limit (3 <C 1, the system behaves like the uncoupled 
system as observed on Figure [TUTa) . 

In the large-/3 limit, the output circuit behaves like an open loop: electric charge transfers in the circuit between 
the two faces of the piezo-electric material are negligible. Because of the piezo-electric coupling, the plate's bending 
induces an electrostatic potential difference in the piezoelectric patches that, in return, creates an additional rigidity 
on the plate [see Eqs. (fTo| - (fl"T)) ]. The system is therefore equivalent to a flexible plate of modified dimensional 
bending rigidity B' = B(l + a 2 ). This behavior is confirmed on Figure [TOT b) where the stability threshold is observed 
to converge at large /3 toward Uq Vl + a 2 , with Uq (A/*) the uncoupled stability threshold. 

Between these two limits /3 <C 1 and j3 3> 1, the minimum velocity required for flutter to develop is observed to be 
significantly increased for low values of M* (heavy flags) while at large M* (typically M* > 50), a destabilization is 
observed and is maximum for j3 = O(l) when the dynamics of the fluid-solid system and the electrical circuit have 
time-scales of the same order. We also observe that for f3 = O(l) the different branches observed on the stability 
threshold for large M* and associated with the instability region of successive high-wave number modes, disappear. 
Beyond the conclusions on the evolution of the stability threshold, these results emphasize the existence of a maximum 
in the coupling effects for values of [3 = O(l). 



D. Conversion efficiency 

1. Definition in the finite-length system 

Total harvested energy estimates are not possible within a purely linear framework. However, the conversion efficiency 
(or normalized harvested energy) defined in Eq. (|23[) provides some important information about the energy transfers 
in the linearly unstable modes of the coupled systems and its ability to produce electrical power from the flapping of 
the flexible plate. The spatial average in Eq. ([23]) is now taken over the entire length of the plate and the conversion 
efficiency r becomes in the non-dimensional notations of this section 



t+7T 



v 2 dt'dx. 
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v 2 dt'dx 
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For a function f(x, t) = Re (F(x)e~ Iwt ) with Cj = &) r + iu>i (f being one of ti, q, w...), 
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and Eq. (|64|) becomes 



Gi +Re 

r{t) = 

Hi +Rc 

with Gi = |G(a;)| 2 da; and G 2 = J* G(x) 2 dx, with 

G{x) = 



UJ 



V{xf 



(64) 



(65) 



(66) 
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FIG. 10: Evolution of the stability threshold when ft is varied for a — 0.5. For clarity, only values of j3 lower (resp. greater) 
than 1 are represented in (a) [resp. (b)]. In (a), the solid black line corresponds to the threshold without any coupling (no 
piezo-electric) . In (b), the solid grey line corresponds to an uncoupled system with modified rigidity B' = B(l + a 2 ). 



and Hi and H2 are defined similarly from 

H(x) = \ ([iuW(x)} 2 + -^W"(x) 2 + \v{xf^j . (67) 

Unlike in the local analysis of Sec. IIII1 and because the different functions are not periodic on the averaging length 
anymore, r is still a function of time, except at the stability threshold where cDj = 0. r(t) is however periodic, and for 
unstable flutter modes, the fluctuations of r around its mean value are weak as the mode's growth rate is generally 
small compared to its frequency. In the following, most of the analysis will be performed at the stability threshold, 
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FIG. 11: Conversion efficiency of the dominant unstable mode for a = 0.5 and ft — 0.25. In the white regions, all the modes 
are stable. The stability threshold obtained for the uncoupled system (a = 0) is shown as a dashed line. 



where r is a constant. In all other cases, its time-average r over a period will be taken as a measure of the conversion 
efficiency and is easily computed from Cj and the mode shape functions V and W. 



We are interested here in the coupled behavior of the fluid-solid and electric systems for a given piezoelectric material 
(fixed a). Figure [Til shows the variations with M* and U* of the conversion efficiency associated with the dominant 
mode for a = 0.5 and /3 = 0.25, above the stability threshold. We observe that this conversion efficiency is increased 
significantly for large M* , corresponding to high fluid-to-solid mass ratios (in water for example) or long systems. The 
conversion efficiency is also observed to be maximum in the regions destabilized by the piezo-electric coupling. This 
is consistent with the conclusions of the local analysis, where a maximum of the conversion efficiency was observed 
for the Negative Energy Wave destabilized by the introduction of piezoelectric damping. 

Such destabilized regions are located in the vicinity of the instability threshold and we will from now on focus on 
this parameter region, as it is also most relevant given the linear character of this study. The map of the conversion 
efficiency for the dominant (and neutrally stable) mode at the threshold is shown on Figure [12] It confirms the 
existence of a maximum of conversion efficiency at high M* for /? ~ 0.1-1. We also observe that the maximum 
efficiency and the corresponding value of j3 are independent of M* for M* > 20-50. This suggests that the dynamics 
leading to the maximum energy conversion are dominated by local effects as will be discussed further in Section [V] 
Figure [12] also confirms that the maximum conversion efficiency is reached when the flexible plate is destabilized by 
the piezo-electric coupling. 

In that region of the parameter space, the frequency of the dominant mode varies between 3 and 5 (Figure fT3]) . 
which is of the same order of magnitude as the value of 1//3 near the maximum efficiency, therefore suggesting a 
resonance-type phenomenon between the frequency of the neutrally stable mode and the characteristic frequency of 
the output electric system. The existence of a maximum for the conversion efficiency for (3u> r ~ 1 is confirmed on 
Figure [T3] Note that the agreement is excellent for M* = 100 while it is not as good for M* = 1 which lies outside 
of the maximum efficiency region. 

The evolution of the maximum efficiency with /3 on Figure Q2] is highly reminiscent of the variations of the dissipated 
energy in the resistive part of a i?C-electrical circuit forced by a generator at a given frequency uj. The role of the 
forcing generator is played here by the fluid-solid system that imposes charge transfers and non-zero potential in the 
piezo-electric material through the deformation of the flexible plate [Eq. (|17[) ]. The forcing frequency is however not 
independent from the output circuit, as shown on Figure [T3] where we observe a significant variation with (3 of the 
frequency of the dominant mode. This modification of the forcing frequency is a result of the piezoelectric feedback of 



2. Energy conversion and resonance 
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FIG. 12: Conversion efficiency r of the dominant mode (neutral) at the stability threshold for a — 0.5. The thick black line 
shows the limits of the regions of the (M* , /3)-plane destabilized in comparison with the uncoupled problem a = 0. 




FIG. 13: Evolution with fi at the stability threshold of the conversion efficiency (thick solid) and frequency (thick dashed) for 
a = 0.5 and M* — 1 (a), M* — 10 (b) and M* = 100 (c). The evolution of the critical velocity U* at the threshold is also 
shown on each plot (thin dash-dotted) and the thin solid line represents 1//3 for reference. The vertical dotted line indicates 
the value of /3 for which the conversion efficiency is maximum. 



the electrical system on the fluttering dynamics. In that regard, the observed efficiency peak differs from a traditional 
resonance and corresponds to a tuning of the fluid-solid frequency to the output system's when the two time-scales 
are of the same order. 

The comparison of the three cases M* = 1, M* = 10 and M* = 100 also confirms that the maximum conversion 
efficiency is significantly larger when destabilization by damping occurs (lower value of the threshold when /3w r ~ 1). 
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FIG. 14: Evolution of the maximum conversion efficiency r max with the coupling coefficient a. r max is defined as the maximum 
of r evaluated at the stability threshold for all values of j3 and M * . The insert plot represents the same quantities in log- log 
format, showing that r max ~ a 2 . 



3. Influence of the piezoelectric coupling 

In Figure [121 we observed that the conversion efficiency reaches a maximum value of 1.7 for large M* and /3u> r ~ 1. 
Figurc[l4]shows that the maximum conversion efficiency at the stability threshold is strongly influenced by the coupling 
coefficient a, with a clear scaling r max ~ a 2 , and this confirms the discussion of Section IlIIl It is not surprising to 
obtain values of r max greater than one as it is not a thermodynamic efficiency but a measure of the harvested energy 
over a period relative to the average stored energy over the same period in the solid and capacitative systems. Note 
also, that unlike /?, U* and M*, that are determined by the tuning of the output system properties, of the plate's 
length and inertia and of the flow conditions, a only depends on the characteristics of the piezo-electric material 
and the plate's rigidity. For a given piezo-electric system a is fixed and Figure [14] provides an upper bound of the 
achievable conversion efficiency. 



V. DISCUSSION 



The fluttering of a flexible plate in an axial flow is an attractive candidate for flow energy harvesting as it can produce 
self-sustained periodic vibrations of the solid body. Understanding the impact of the harvesting-induced damping on 
the fluid-solid dynamics is an essential element to assess the potential performance of such systems. Thus the present 
paper investigated the linear dynamics of the fluttering flexible plate fully coupled to a simple dissipative electrical 
circuit through piezoelectric layers, by studying the coupling effect on the local and global instabilities as well as on 
the energy conversion efficiency. 

The ability to quantify the total harvested energy is intrinsically limited by the linear framework considered, as 
saturated flapping amplitudes are not computed. However, the linear analysis is an important step to obtain some 
insight on the modification of the system properties, in particular its stability threshold, by the introduction of a 
realistic coupling to an output electrical network. Hence, it was shown that destabilization by damping can occur 
in such systems. The local analysis showed that the destabilized modes correspond to the so-called Negative Energy 
Waves of the un-damped system. This was confirmed in the global analysis by the destabilization of plates with high 
mass ratio where local effects arc expected to dominate. In particular, we observed on Figure [121 that above a certain 
value of M* of the order of 100, the dependence with M* of the system properties at the stability threshold becomes 
negligible. 
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The parameter M* can be seen as a relative measure of the length of the system to the local characteristic lcngthscale. 
For M* 3> 1, one expects the stability and properties of the dominant modes to be driven by local phenomena rather 
than global ones. For given a and ft, the critical velocity threshold remains finite which is consistent with the 
observation that the system is unstable locally for all V* = U*M*. The presence of a maximum in conversion 
efficiency for ftu r ~ 1 can be identified to the maximum efficiency observed in local analysis for ^uj r ~ 1 as ft = j/M* 
and uj r = Cj t M* . 

From a practical point of view, the present results suggest that the energy conversion is more efficient when the mass 
ratio M* is large. This is achieved when the plate's length is long, but also when the fluid inertia is large compared 
to the solid's, as in water flows for example. These results also emphasize the importance of the output circuit in 
the energy transfers: a careful tuning of the circuit characteristic time-scale to that of the fluid-solid oscillations 
significantly increases the conversion efficiency from the solid to the electric system. 

Using both local and global analyses, the maximum energy conversion efficiency r max was found to scale as a 2 . 
This scaling is expected in the limit of weak coupling: if one neglects the piezo-electric coupling term in Eq. (f24)) or 
Eq. (|52]l. the fluid-system can be considered as only weakly modified and act as a constant forcing on the electric 
system. The forcing in Eq. (|25|) and Eq. (|53| scales linearly with a and so is expected to scale the charge density. 
The harvested energy, a quadratic function of the charge density, is therefore expected to scale as a 2 . It is however 
surprising to observe that this scaling remains valid in a greater range of coupling. This result enlightens the essential 
role of the coupling coefficient a in energy harvesting applications. This coupling coefficient can be obtained from 
Eqs. ^ and (f2"T)) by separating the properties of the piezoelectric materials from the geometric effects of the relative 
thickness of the plate and piezoelectric patches: 

\ l ~ V l cs( h * /- 1 




hp E p (l - v%) 



where §f is a non-dimensional function of the thickness and bending stiffness ratios. To increase the piezoelectric 
coupling it is therefore important to maximize e^i or minimize the permittivity s. It must be noted that changing 
e also impacts the capacity of the piezoelectric material, which in turn modifies the choice of the optimal energy 
harvesting electrical circuit through the optimization of the parameter ft. Finally, once the materials are chosen, 
an optimal thickness ratio ho/h p can also be obtained. To get an idea of the value of a in practical applications, 
two cases may be considered that are representative of typical values found in the literature: the former consists 
in a mylar plate (ho = lOO^im, Eq = 4GPa) with two PVDF piezoelectric layers (h p = 40/^m, E p = 2.5GPa, 
e 31 = 0.023C/m 2 , e r = 11.5), leading to an approximate value of a ~ 0.03. The latter consists in a steel plate 
(ho = 300/im, Eo = 200GPa) with two PZT piezoelectric layers (h p = 300/im, E p = 60GPa, e 31 = 10C/m 2 , 
e r = 2000) and in that case, the coupling coefficient is a ~ 0.3. Careful design of the system is expected to further 
increase a, therefore the typical value of a = 0.5 chosen in the present paper is realistic but corresponds to an upper 
bound estimate of the system's performance. 

Natural extensions of this work include the study of passive resonant [32j or active circuits [43[ to determine 
potential improvements in the energy transfer from the electrical design. In this paper, a continuous distribution of 
piezoclcctrics has been considered and studying the impact of the finite-length of piezoelectric patches would provide 
important information on the actual design of energy harvesters. Finally, as emphasized throughout the present 
paper, further investigation of the energy harvesting potential must include the representation of non-linear effects 
in the fluid and solid dynamics, to obtain the amplitude of the self-sustained oscillations of the system. Both local 



44j and global [201 ] analyses of the nonlinear regime have been performed for compliant panels or plates placed in an 
axial flow. The- non-linear behavior of piezoelectric materials in the energy harvesting context has been investigated 
by [45| . The study of the saturated regime in the fluttering dynamics of the fluid-structure-electrical fully-coupled 
system will provide a quantitative assessment of the actual harvested power and is the focus of subsequent work to 
this linear analysis. 
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